\documentclass[%
    %draftclsnofoot,
    %submission,
    %compressed,
    %final,
    %
    %technote,
    %internal,
    %submitted,
    %inpress,
    %reprint,
    %
    %titlepage,
    %notitlepage,
    %anonymous,
    %narroweqnarray,
    %inline,
    twoside,
    %    invited,
    ]{IEEEtran}
%\bibliographystyle{IEEEtran}


% *** MATH PACKAGES ***
%
\usepackage[fleqn]{amsmath}
\usepackage{amsfonts}
\usepackage{graphicx}
%\usepackage{tikz}
 %\graphicspath{{./eps/}}
%\usepackage{subfigure}

\usepackage{soul}
\usepackage{threeparttable}


%
%% correct bad hyphenation here
%\hyphenation{op-tical net-works semi-conduc-tor}

\begin{document}
%
% paper title
% can use linebreaks \\ within to get better formatting as desired
\title{The Fokker-Planck Equation for Power System Stability Probability Density Function Evolution}


\author{Keyou~Wang,~\IEEEmembership{Member,~IEEE,}
        Mariesa~L.~Crow,~\IEEEmembership{Fellow,~IEEE}
% <-this % stops a space

%\thanks{The authors are with the Department of Electrical and Computer Engineering, Missouri University of Science and Technology, Rolla, MO 65401 USA.}% <-this % stops a space

%%\thanks{Manuscript received April 19, 2005; revised January 11, 2007.}
}

% make the title area?
\maketitle


\begin{abstract}
%\boldmath
This paper presents an analysis of the evolution of the probability density function of the dynamic trajectories of a single machine infinite bus power system.  The probability density function can be used to determine the impact of random (stochastic) load perturbations on system stability. The evolution of the state probability density function over time leads to several interesting observations regarding stability regions as a function of damping parameter.
The Fokker-Planck equation (FPE) is used to describe the evolution of the probability density of the states. The FPE is solved numerically using standard partial differential equation (PDE) solvers (such as finite difference method).  Based on the results, the qualitative changes of the stationary density produce peak-like, ridge-like and other complicated shapes. Lastly, the numerical FPE solution combined with single machine infinite bus (SMIB) equivalent techniques lay the framework extended to the multimachine system.

%In addition, the numerical results indicate that under certain conditions additive random noise may actually improve the stability of system.

%
%
%development of the exact stationary solution of the Fokker-Planck equation (FPE) for the nonlinear multi-machine power system under white noise excitation. The proposed solution of FPE describes the stationary probability density function (PDF) of the states after the system reaches statistical stationarity. The analytical result shows that the stationary PDF of the stochastic power system is a function of the Hamiltonian. The proposed solution is illustrated with a single machine power system and multiple machine power system. The proposed exact stationary solution is also validated through a comparison with several numerical solutions of FPE including the finite difference method and a Monte Carlo approach.
\end{abstract}

% Note that keywords are not normally used for peerreview papers.
\begin{IEEEkeywords}
Fokker-Planck equation, Stationary stochastic processes, Probability density function, Stochastic differential equations, Finite difference methods, Power system stability
\end{IEEEkeywords}


% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
%\IEEEpeerreviewmaketitle

\newtheorem{theorem}{Theorem}
\newtheorem{remark}{Remark}
\newtheorem{definition}{Definition}
\newtheorem{proposition}{Proposition}
\newtheorem{lemma}{Lemma}
\newtheorem{assumption}{Assumption}


\section{Introduction}
\PARstart{B}{y} nature, a power system is continually experiencing stochastic disturbances, whether as  switching events or through generation and load level fluctuations \cite{fp:Kundur94}. The switching event disturbances include generator outages, short-circuits caused by lightning or other fault conditions, sudden generator outages, large load changes, or a combination of such events.  Many event disturbances lead to a change in the topology of the power system between pre-fault and post-fault conditions.  A disturbance which causes a change in topology may suffer an immediate change in stability since the post-fault equilibrium may have different stability properties than the pre-fault system.  On the other hand, small random disturbances of power generation and load may gradually change the region of stability over time. In today's utility environment, the variability of electrical network loading has received increased attention in recent years due to the expansion of  renewable energy resources such as wind turbines and solar panels, and the likelihood of wide-spread adoption of plug-in electric vehicles (PEVs).  This tandem effect of renewable resources and PEVs may create uncertainties of such significant magnitude they impact the operation of the power system and make transient stability assessment difficult.


%\PARstart{P}{ower} system transient stability is always the one of the kernel concerns for the power industry. In past decades, power system stability assessment has been accomplished by the time-domain simulation or through the energy-function based direct method \cite{fp:Kundur94}. Traditional approaches have both been developed based on deterministic models and have served the power industry reasonably due to light uncertainties of the loading and generation. However, in today's environment, the expansion of renewable energy resources create uncertainties of such significant magnitude they impact the operation of the power system and cannot be neglected. There is a need to account for the probabilistic nature of system conditions and events and to quantify and manage the risk.
%

Previous transient stability stochastic studies have addressed uncertainty related to disturbances and operating conditions through a combination of deterministic time-domain simulation techniques with stochastic analyses \cite{fp:Billintion81}-\cite{fp:Hockenberry04}.  The framework of these studies has been to numerically integrate the differential equations and then iteratively repeat the integration procedure to yield the statistical information through Monte Carlo methods \cite{fp:Billintion81}-\cite{fp:Wu83}, the probabilistic collocation method \cite{fp:Hockenberry04}, or the trajectory sensitivity method \cite{fp:Hiskens06} to obtain the first and second moment values (mean and variance) of the states.  In these approaches however, the evolution of the probability density function (PDF) over time has not been studied.

%

Although it is well-known that the numerical solution of stochastic differential equations (SDEs) is a useful tool to study the trajectory of stochastic power system, little work has focused on the evolution of the probability density function during the transient period. In this paper, we shift our perspective from the sampling trajectories of SDEs to the probability density distribution. The Fokker-Planck equation (FPE) is a parabolic partial differential equation (PDE) which can be derived from the corresponding SDE and describes the evolution of probability density of the state variables in the stochastic system \cite{fp:Risken89}. Our intention is to apply the Fokker-Planck equation to the study of the probability density function in the context of stochastic power system stability.

%
In this paper, we start our analysis by considering the swing equation for a single machine connected to an infinite bus (SMIB) system perturbed with stochastic excitation. The dynamics of the deterministic classical SMIB system has been extensively studied \cite{fp:Kundur94}, \cite{fp:Paganini99}-\cite{fp:Chen05}.  The SMIB model has an inherent relation with many physical systems such as the forced pendulum \cite{fp:Andronov66} or the Josephson junction circuit \cite{fp:Levi78}, which both exhibit two attracting equilibria including a single stable equilibrium point and a stable limit cycle.  We will build upon these deterministic analyses to develop predictions for the behavior of stochastically perturbed power systems.   Specifically, the primary contributions of this paper are to
\begin{enumerate}
\item analyze the stability properties of the stochastic SMIB,
\item develop and solve the corresponding Fokker-Planck equation to predict the transition of the system states,
\item quantify the time-dependent probability of stability, and
%\item draw correlations and conclusions regarding the behavior and stability of the deterministic versus stochastic power system.
\item combine equivalent SMIB techniques to lay the framework of multimachine stochastic systems.
\end{enumerate}

\hl{Especially extended by the SMIB equivalent techniques, the proposed methodology therefore have potential to provide the system operators with quantitative measure of probability of transient stability in practical large-scale multi-machine systems under stochastic perturbations.}
%
%In this paper, we review the numerical solution of stochastic differential equation and use of monte carlo approach to approximate the probability density. Then we propose the Fokker-Planck equation (FPE) to directly describe the probability density of the stochastic SMIB model with white noise excitation. An implicit finite-difference scheme is proposed to solve the corresponding FPE numerically. Numerical results are presented to provide insights to understand the evolution of probability density upon parameters.

\section{Stochastic SMIB System}

To provide a framework for developing the probability density function, a stochastic SMIB system is initially considered.  The normalized SMIB model with additive white noise perturbations is modeled as a set of Ito stochastic differential equations (SDE) \cite{fp:Oksendal07},
  \begin{eqnarray}
    \!\!\!d \delta&=& \omega d t \label{eq:SMIB SDE 1}\\
     \!\!\!d \omega &=& \left( \bar{P}_m-\sin \delta-\bar{D}\omega \right)d t + \sigma dW(t) \label{eq:SMIB SDE 2}
  \end{eqnarray}
where $\delta$ is  the relative generator angle to the infinite bus, $\omega$ is the generator velocity (relative to synchronous speed). The state vector is  $x=[x_1, x_2]=[\delta, \omega]$. The constants $\bar{P}_m$ and $\bar{D}$ are the normalized mechanical power and normalized damping coefficients, respectively. The function $W(t)$ is a Wiener process corresponding to the Gaussian white noise excitation dW(t). The coefficient $\sigma$ denotes the magnitude of the additive white noise. This noise can be considered as the stochastic disturbance of the bus load power. \hl{Typically, SDEs incorporate white noise which can be thought of as the derivative of Brownian motion (or the Wiener process); however, it should be mentioned that other types of random fluctuations are possible.}
%i.e. $\bar{P}_m=\bar{P}^0_m+\sigma \frac{dW}{dt}$.


\hl{Gaussian white noise is natural choice for modeling many random variables which will be normally distributed and, such as random load level in power system}\cite{fp:Odun-Ayo12},\cite{fp:Dong12}. \hl{ However, it should be noted that other more complex random fluctuations can be incorporated.} \hl{Poisson process is suitable to model the discrete random events, which are able to represent the system faults}\cite{fp:Dong12},\hl{ or sudden change of output of load/generation.} In \cite{fp:Wang12}, \hl {the Rayleigh process can be constructed by an one-dimensional SDE which has the stationary PDF of Rayleigh distribution. And Rayleigh distribution is a particular Weibull distribution that is well-known to model wind energy.}

%The novelty of our study is that

\section{Numerical Solution of Stochastic Differential Equation }

It is computationally difficult to study the time-dependent analytical solution of (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}) over the series of Monte Carlo trials necessary to produce meaningful results. The numerical solution of SDEs requires the time-domain numerical integration of (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}) via weak or strong methods \cite{fp:Kloeden99}. The primary difference between the numerical solution of the system of ODEs and SDEs is  the additional inclusion of the standard Wiener process over the simulation interval $[0, T]$.  The Wiener process $W(t) (t\geq 0)$ is a stochastic process that satisfies the following conditions \cite{fp:Higham01}.
\begin{enumerate}
\item $W(t)$is continuous and $W(0)=0$ (with probability 1);
\item $\forall t\geq 0, W(t)\sim \sqrt{t}\mathcal{N}(0,1)$
\item $W(t+\delta t)-W(t) \sim \sqrt{\delta t} \mathcal{N}(0,1)$  is
independent of the history of the process up to time $t$
\end{enumerate}
where $\mathcal{N}(0,1)$ denotes a normally distributed random variable with zero mean and unit variance.

%\subsection{Euler Maruyama Scheme}

Writing (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}) in  vector form gives:
  \begin{equation}
dx = f(x) d t + g(x) dW(t) \label{eq:SMIB SDE 3}
  \end{equation}
where
\begin{eqnarray}
   f(x)&=& \left[ \begin{array}{c}
     f_1(x)\\
     f_2(x)
   \end{array}\right]=\left[ \begin{array}{c}
     x_2\\
     \bar{P}_m-\sin x_1-\bar{D}x_2
   \end{array}\right] \\
   g(x)&=& \left[ \begin{array}{c}
   0\\
   \sigma
   \end{array}\right]\\
   dW(t)&=&W(t+\delta t)-W(t)
\end{eqnarray}

To numerically integrate (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}) over the interval $[0, T]$, the interval is discretized with $0=\tau_0 < \tau_1<\cdots <\tau_n<\cdots<\tau_N=T$. The equidistant case has step size $\Delta t=\frac{T}{N}$, and $\tau_n=n \Delta t$, where $x_n$ is the discrete approximation of the exact solution $x(\tau_n)$. Then the Euler-Maruyama (EM) integration scheme takes the form \cite{fp:Higham01}:
  \begin{equation}
      x_{n+1} = x_n+ f(x_n) \Delta t + g(x_n) (W(\tau_{n+1})-W(\tau_{n})) \label{eq:EM Scheme}
  \end{equation}
for $n=0,1,\ldots,N-1$ with initial value $x_0=x(0)$. The values $W(\tau_{n+1})$,$W(\tau_{n})$ are points on the discretized Brownian path which are generated for the specific integration interval. The Euler-Maruyama method is used because it is the simplest strong Taylor appropriation, containing only the time and Wiener integrals of multiplicity one from the Ito-Taylor expansion. The EM scheme usually attains the order of strong convergence of $\frac{1}{2}$, whereas the underlying deterministic Euler method converges with classical order 1. Other SDE methods can be applied, such as the Milstein's higher order scheme, to raise the strong convergence order. Moreover, the implicit stochastic schemes are also desirable for the sake of numerical stability as well as deterministic schemes.

%\begin{figure}
%\centering
%  \hspace*{\fill}{\includegraphics[width=3.5in]{SMIB.eps}}
%\caption{Deterministic SMIB system}
%\label{fig:SMIB ODE}
%\end{figure}
%
%
%To illustrate the numerical solution of the stochastic differential equations, the method is applied to solve (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}). As a benchmark, the deterministic system is also solved numerically to provide  a benchmark for the stochastic results (the deterministic ordinary differential equations are obtained by setting $\sigma=0$ in (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2})).  The initial conditions are selected to lie in the attraction region of the stable equilibrium point (SEP) of the deterministic system, but close to the stability boundary. In the deterministic case, the trajectory is expected to converge to the SEP. The resulting generator angle and frequency are shown in Fig. \ref{fig:SMIB ODE}.

\begin{figure}
\centering
\includegraphics[width=2.6in]{sde_delta.eps}
\caption{Stochastic SMIB system results (20 runs) - generator angle}
\label{fig:SMIB SDE delta}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=2.6in]{sde_w.eps}
\caption{Stochastic SMIB system results (20 runs) - generator frequency}
\label{fig:SMIB SDE w}
\end{figure}


To illustrate the numerical solution of the stochastic differential equations, the method is applied to solve (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}). The initial conditions are selected to lie in the attraction region of the stable equilibrium point (SEP) of the deterministic system, but close to the stability boundary. In the deterministic case, the trajectory is expected to converge to the SEP. However, in Figs. \ref{fig:SMIB SDE delta}-\ref{fig:SMIB SDE w}, twenty stochastic generator angle and frequency trajectories are shown for different noise sets (i.e. 20 Monte Carlo trials). The majority of the trajectories closely mimic the deterministic dynamic trajectories.  However, two of the twenty trajectories exhibit unexpected behaviors and become unstable.  A mathematical analysis of the two errant trajectories indicate that actually these two trajectories form a ``rotating orbit.''  If the angle were plotted on a rotating reference frame (repeating every $2\pi$ radians), the two irregular phase portraits shown in Fig.  \ref{fig:sde_phase2} would correspond to a limit cycle \hl{ on a cylinder}. \hl{ The portrait is always $2\pi$-periodic in $\delta$ and thus can be folded into a cylinder, as shown in Fig.} \ref{fig:sde_phase_cylinder}. \hl{The dynamics of SIMB system is another analog of the driven pendulum or the Josephson junction circuit. These physical systems all exhibit intriguing hysteresis effects, thanks to the coexistence of a stable equilibrium point and a stable limit cycle.}\cite{fp:Strogatz01}.


\begin{figure}
\centering
\includegraphics[width=2.6in]{sde_phase2.eps}
\caption{Phase portrait showing convergence to SEP and two rotating solutions}
\label{fig:sde_phase2}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=2.6in]{cylinder.eps}
\caption{Phase portrait on a cylinder}
\label{fig:sde_phase_cylinder}
\end{figure}

%\subsection{Monte Carlo Simulation}

Whereas each individual run provides the behavior of a single trajectory, a Monte Carlo approach with multiple trials is  needed to yield statistical information regarding the stochastic behavior of the system.  Note that there are two vertical lines labeled $p(x,t)$ drawn at 5 and 15 seconds in Figs. \ref{fig:SMIB SDE delta}-\ref{fig:SMIB SDE w}.  The cross section of the trajectories denoted by the vertical lines approximates the probability density function (PDF) of the state $x$ (where $x$ is $\delta$ or $\omega$) at  time $t$.  The PDF $p(x,t)$ describes the probable distribution of each state at time $t$.  As the number of trials $M$ increase,
the Monte Carlo method will provide an increasingly accurate approximate of the true probability distribution function. However, since $M$ is a finite number, only an approximate PDF can be obtained.  The use of histograms provides a straightforward approach to estimating (and visualizing) the joint PDF.   Figs. \ref{fig:SMIB SDE_HIST_t1} and \ref{fig:SMIB SDE_HIST_t2} show the approximated probability density at the time $t=5$ s and $t=25$ s respectively for 1000 runs. Consider first Fig. \ref{fig:SMIB SDE_HIST_t2} which shows the histogram of state values at $t=25$ seconds.  This histogram shows two distinct areas of cumulative states:  the large number of states clustered around the SEP (near $(0.5,0)$) and the arc of values encircling the SEP.  This arc (or ridge) corresponds to the set of states that lie on the limit cycle encircling the SEP.  The same distribution is apparent in Fig. \ref{fig:SMIB SDE_HIST_t1} except that at $t=5$ seconds, the states have not yet settled near the SEP.  A series of histogram ``snapshots'' at each second would clearly illustrate the transition of the PDF from $t=0$ to steady-state and how the states are attracted by either the SEP or the limit cycle.

%More sophisticated schemes use kernel-based density estimation in which each sample point is replaced by a small probability density function (such as a Gaussian distribution)\cite{fp:Molina-Garcia11}.

%Each of these is called a kernel. Then the each kernel can contribute together to estimate the overall .

Although the histogram approach can provide a reasonable approximation of the evolution of the PDF, Monte Carlo trials has the disadvantage of being computationally burdensome. Furthermore, due to the discrete nature of the histograms, subtle dynamic features may be overlooked or discarded as ``numerical noise'' resulting from the numerical integration process. Therefore, it is desirable to develop a method that can accurately predict the transient evolution of the PDF and identify salient dynamic features.  The Fokker-Planck equation provides one such robust technique.

%
\begin{figure}
\centering
\includegraphics[width=2.6in]{sde_hist_t1.eps}
\caption{Histogram of PDF at $t=5$ sec }
\label{fig:SMIB SDE_HIST_t1}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=2.6in]{sde_hist_t2.eps}
\caption{Histogram of PDF at $t=25$ sec }
\label{fig:SMIB SDE_HIST_t2}
\end{figure}
%
%\begin{figure}
%\centering
%\includegraphics[width=1.7in]{sde_hist_t1.eps}
%\includegraphics[width=1.7in]{sde_hist_t2.eps}
%\caption{Histogram of PDF at $t=5$ sec }
%\label{fig:SMIB SDE_HIST_t1}
%\end{figure}





\section{Fokker-Planck equation for stochastic SMIB model}


In this section,  a second approach to obtaining the relevant probability density functions is proposed.  It is well-known that if a diffusion process (e.g. a continuous Markov process) is the result of solving the system of stochastic differential equations, then the transitional probability density function $p(x,t|x_0,t_0)$  satisfies the {\it Fokker-Planck} equation (FPE) \cite{fp:Risken89}-\cite{fp:Lin04}.  In this section,
the Fokker-Planck equation is derived from the stochastic differential equations of (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}) to describe the time evolution of the probability density $p(x,t|x_0,t_0)$ for the random variable $x(t)$:

\begin{eqnarray} \label{eq:FPK G}
    \frac{\partial p}{\partial t}=-\nabla \cdot G(x,p)=-\sum\limits_{i=1}^{n}\frac{\partial }{\partial x} G_i
\end{eqnarray}
where {\it probability flows} $G=(G_1,\ldots, G_n)$ are defined by
  \begin{eqnarray} \label{eq:FPK G Def}
    G_i=a_i(x)p(x)-\frac{1}{2}\sum\limits_{j=1}^{n}\frac{\partial}{\partial x_j}\left[ b_{ij}(x) p(x)\right]
  \end{eqnarray}
where $a(x,t)=[a_i]_{n \times 1}$ and $b(x,t)=[b_{ij}]_{n\times n}$ are the {\it drift} coefficients and {\it diffusion} coefficients at time $t$ and position $x$.


In the SMIB case, $a(x,t)$ and $b(x,t)$ can be derived from the Ito stochastic differential equation of SMIB system (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}):
\begin{eqnarray}
   a(x)= f(x)=\left[ \begin{array}{c}
     x_2\\
     \bar{P}_m-\sin x_1-\bar{D}x_2
   \end{array}\right]
\end{eqnarray}
and
\begin{eqnarray}\label{eq:FPK diffu}
 b(x)=g(x) g^T(x)=\left[ \begin{array}{cc}
     0&0\\
     0&\sigma^2
   \end{array}\right]
\end{eqnarray}

From (\ref{eq:FPK G})-(\ref{eq:FPK diffu}), the corresponding Fokker-Planck equation of (\ref{eq:SMIB SDE 1})-(\ref{eq:SMIB SDE 2}) is given by:
%\begin{eqnarray} \label{eq:SMIB FP}
% \frac{\partial p}{\partial t}=-x_2\frac{\partial p}{\partial x_1}-\left( \bar{P}_m-\sin x_1-\bar{D}x_2\right)\frac{\partial p}{\partial x_2}+\frac{\sigma^2}{2} \frac{\partial^2p}{\partial x_2^2}
%\end{eqnarray}
 \begin{eqnarray} \label{eq:SMIB FPK}
 \frac{\partial p(x,t)}{\partial t}&=&-\frac{\partial f_1(x)p(x,t)}{\partial x_1}-\frac{\partial f_2(x)p(x,t)}{\partial x_2} \notag \\
 &&+\frac{\sigma^2}{2} \frac{\partial^2p(x,t)}{\partial x_2^2}
\end{eqnarray}
where $p(x,t)$ is  the time-dependent probability density function. The FPE (\ref{eq:SMIB FPK}) is a two-dimensional parabolic partial differential equation with specified initial conditions and boundary conditions.

The \emph{initial condition} of (\ref{eq:SMIB FPK}) is given by
\begin{equation}\label{eq:FPK Ini}
 p(x_0|x,t)=\prod\limits_{i=1}^{n} \tilde{\delta}(x_i-x_{i0}) \quad \text{as} \quad t\rightarrow 0
\end{equation}
where $x_0$ is the initial state and $\tilde{\delta}$ is the Dirac delta function.

Solutions of the nonlinear FPE (\ref{eq:SMIB FPK}) depend significantly  on the boundary conditions \cite{fp:Lin04}. Accordingly, two sets of boundary conditions are imposed on the SMIB problem:

\begin{figure}[b]
\centering

%\def\AA{-2.2}
%\def\BB{2.2}
%\def\xmin{-1.8}
%\def\xmax{1.8}
%\begin{tikzpicture}[domain=\AA:\BB]
%\draw[color=gray] plot[smooth] (\x, {exp(-1*\x*\x)});
%\filldraw [fill=gray!30] plot[id=f1,domain=\xmin:\xmax] (\x, {exp(-1*\x*\x)})
%            -- (\xmax,0) -- (\xmin,0) -- cycle ;
%\draw[dashed] (\xmin,1.2) -- (\xmin,0) node[below] {$x_{2min}$};
%\draw[dashed] (\xmax,1.2) -- (\xmax,0) node[below] {$x_{2max}$};
%\draw[->] (\AA,0) -- (\BB,0) node[right] {$x_2$};
%\draw[->] (0,-0.2) -- (0,1.5) node[above] {$p(x_1,x_2,t)$};
%\end{tikzpicture}
\includegraphics[width=2.6in]{fig6.eps}
\caption{Approximation of natural boundary conditions in the direction of $x_2$}
\label{fig:x2 boundary}
\end{figure}

\begin{enumerate}
  \item \emph{Natural boundary conditions:}
  \begin{eqnarray}
      p(x_1,+\infty,t) &=& 0\\
      p(x_1,-\infty,t) &=& 0
  \end{eqnarray}
These boundary conditions require that the probability density should vanish at $x_2 \rightarrow -\infty$ and $x_2 \rightarrow +\infty$. In the numerical solutions of the corresponding FPE (\ref{eq:SMIB FPK}), two finite value are taken to approximate the boundary for $x_2$, as illustrated in Fig. \ref{fig:x2 boundary}. Therefore the natural boundary conditions can be replaced by the { \it absorbing  boundary conditions}:
  \begin{eqnarray}
      p(x_1,x_{2min},t) &=& 0\\
      p(x_1,x_{2max},t) &=& 0
  \end{eqnarray}
where the $x_{2max}$ and $x_{2min}$ are finite and define the boundary of computational area along the $x_2$ direction.


  \item \emph{Periodic  boundary conditions:}
  \begin{eqnarray}
         p(x_1,x_2,t) &=& p(x_1+2\pi,x_2,t)
  \end{eqnarray}
Since $f(x)$ is periodic with $x_1$ ($\delta$) in $2\pi$, the  state space $(x_1, x_2)$ can be viewed either in planar $R^2$ or cylindrical $S^1 \times R$. For simplicity, the cylindrical state-space $S^1 \times R$ is assumed for the solution of the FPE (\ref{eq:SMIB FPK}) to capture the near periodic behavior of the limit cycle.
\end{enumerate}

A normalization condition for the probability density function  must also be satisfied:
\begin{eqnarray}
  \int p(x,t) dx &=& 1
\end{eqnarray}

For the most part, the steady-state solution (also called the ``stationary solution'') of the FPE is of greatest interest:  i.e. the behavior of the system dynamics as the states approach steady state. Note the nonlinear FPE may have no stationary solution. If the stationary solution exists as $t \rightarrow \infty$, then $p_\infty(x|x_0)$ is the stationary probability density function. In other words, $p_\infty(x|x_0)=\lim_{t\longrightarrow \infty} p(x,t|x_0)$. Under this steady-state condition, (\ref{eq:SMIB FPK}) becomes the stationary FPE where:
\begin{eqnarray} \label{eq:FPK G Stationary}
    \sum\limits_{i=1}^{n}\frac{\partial }{\partial x} G_i =0
\end{eqnarray}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time-dependent pdf, Pm=0.5, D=0.409, sigma=1
\begin{figure*}
\centering
%\begin{minipage}[b]{0.32\textwidth}
%\centering
%\subfigure[t=0]{
\label{fig:td pdf:a}
\includegraphics[width=1.8in]{tdpdf-surf-t1.eps}
%\end{minipage}%
%\begin{minipage}[b]{0.32\textwidth}
%\centering
%\subfigure[t=5]{
\label{fig:td pdf:b}
\includegraphics[width=1.8in]{tdpdf-surf-t2.eps}
%\end{minipage}
%\begin{minipage}[b]{0.32\textwidth}
%\centering
%\subfigure[t=25]{
\label{fig:td pdf:c}
\includegraphics[width=1.8in]{tdpdf-surf-t3.eps}
%\end{minipage}
%
%
%\begin{minipage}[b]{0.32\textwidth}
%\centering
%\subfigure[t=0]{
\label{fig:td pdf:e}
\includegraphics[width=1.8in, height=1.35in]{tdpdf-cont-t1.eps}
%\end{minipage}
%\begin{minipage}[b]{0.32\textwidth}
%\centering
%\subfigure[t=5]{
\label{fig:td pdf:f}
\includegraphics[width=1.8in, height=1.35in]{tdpdf-cont-t2.eps}
%\end{minipage}
%\begin{minipage}[b]{0.32\textwidth}
%\centering
%\subfigure[t=25]{
\label{fig:td pdf:g}
\includegraphics[width=1.8in, height=1.35in]{tdpdf-cont-t3.eps}

\begin{scriptsize}
\begin{tabular}{ccc}
\hspace{0.15in} $t=0$ seconds & \hspace{1in} $t=5$ seconds \hspace{1.15in} &$t=25$ seconds
\end{tabular}
\end{scriptsize}
%\end{minipage}

\caption{Time-dependent probability density: (top) 3D $p(x_1,x_2,t)$, (bottom) 2D contours}
\label{fig:td pdf}
\end{figure*}

\section{Solution of the Fokker-Planck Equation}
\subsection{Time-dependent density functions}

The time-dependent density function $p(x,t|x_0)$ describes the evolution of the transient density function with respect to time $t$. It can be calculated numerically (using a finite difference method) \cite{fp:LeVeque07}. Fig. \ref{fig:td pdf} shows the numerical results of the time-dependent density function three-dimensional contours and the projection of the contours on a planar surface for three different snapshots at $t=0, 5$ and 25 seconds.
%with the following parameters $\bar{P}_m=0.5$, $\bar{D}=0.3$ and $\sigma=0.1$. Here we chose the sequence $t=0,5,25$ for the snapshots.

At $t=0$, the set of initial values must satisfy (\ref{eq:FPK Ini}). In Fig. \ref{fig:td pdf}, the initial state $x_0=[0, 2]$ is chosen. The initial condition is easily seen from the $t=0$ contours.

At $t=5$ seconds, the time-dependent probability density functions are undergoing transient transformations. The majority of the probability flows are approaching the SEP, but some of them have left the main flow and are approaching the limit cycle ridge. After a large time interval, such as $t=25$ seconds in this case, the probability density function converges to its stationary solution. There are two salient features in the stationary solution: one is a peak surrounding the SEP, whereas the other is a ridge defining the limit cycle.

Since the stability of the system is predicted by the stationary probability density function, the remainder of the paper will focus on the analysis of the stationary PDF.

\begin{figure}
\centering
\includegraphics[width=2.2in]{spdf-surf-D1.eps}
\caption{Stationary probability density of a highly damped system}
\label{fig:spdf varyD1}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=2.2in]{spdf-surf-D2.eps}
\caption{Stationary probability density of a lightly damped system}
\label{fig:spdf varyD2}
\end{figure}


 \begin{figure}[h]
\centering
\includegraphics[width=2.6in]{logSindex.eps}
\caption{The plot of $\ln(S_{index})$ vs $\bar{D}$}
\label{fig:spdf logS}
\end{figure}

\subsection{Stationary probability density functions}



The stationary probability density function $p_\infty(x|x_0)$ refers to the solution of the FPE after it has reached steady-state. In a highly damped (e.g. highly dissipative) system, the stationary PDF will be centered about the SEP.  However in systems that are lightly damped, the stationary PDF will have a much greater proportion of probability flow along the limit cycle ridge.  Figs. \ref{fig:spdf varyD1} and \ref{fig:spdf varyD2} show examples of the two extremes.  In Fig. \ref{fig:spdf varyD1}, the probability density function is tightly clustered in a peak around the SEP.  This indicates that nearly all trajectories will approach the SEP with high probability.  Fig. \ref{fig:spdf varyD2} illustrates the other extreme of a very lightly damped system in which nearly all of the trajectories are captured by the limit cycle.

In the general case, the stationary probability density function will have a combination of peak and ridge characteristics which depend on the system parameters.  For example, if $\bar{P}_m$ is fixed and $\bar{D}$ is varied, the height of the peak \hl{increases} as $\bar{D}$ increases, whereas the height of the ridge \hl{decreases}. The peak corresponds to the stable case in the context of the SMIB.  The ridge arises from dynamics that are typically considered to be unstable, since the limit cycle only exists for a continually increasing generator angle.



\begin{figure*}
\centering

\centering

\label{fig:SMIB_solutions:a}
\includegraphics[width=1.8in]{pp3.eps}
\includegraphics[width=1.8in]{pp2.eps}
\includegraphics[width=1.8in]{pp1.eps}
\begin{scriptsize}
\begin{tabular}{ccc}
(a) $\bar{D}>\bar{D}_c$, & \hspace{1.0in} (b) $\bar{D}=\bar{D}_c$,\hspace{1.1in} & (c) $\bar{D}<\bar{D}_c$,
\end{tabular}
\end{scriptsize}
\caption{Phase portraits of SMIB system.}
\label{fig:SMIB_solutions}
\end{figure*}


A quantitative measure $S_{index}$ is therefore proposed to relate the stationary PDF to the damping coefficient:
 \begin{eqnarray}
 % \nonumber to remove numbering (before each equation)
   S_{index} &=& \frac{p^{p}_{max}}{p^{r}_{max}}
 \end{eqnarray}
where $p^{p}_{max}$ and $p^{r}_{max}$ are the local maximum values of stationary probability density function  $p_\infty(x)$ near the peak and the ridge, respectively. The value $p^{p}_{max}$ corresponds to the SEP of the deterministic system. The exact location of $p^{r}_{max}$ on the ridge is difficult to determine exactly because of its nearness to the limit cycle, but can be estimated from the numerical solution.   Fig. \ref{fig:spdf logS} shows a series of data points that were estimated at different damping levels.  Note that their relationship can be linearly approximated (with $\ln(S_{index})$) for fixed parameters $\bar{P}_m=0.5$ and $\sigma=0.1$. This allows general conclusions to be drawn regarding the damping level and the resulting dynamic system behavior.  The deterministic critical damping level is calculated to be $\bar{D}_c=0.409$.  In the deterministic system, if the damping parameter is greater than $D_c$, then the system will exhibit only a single SEP and no limit cycles.  However, if the damping is decreased below $D_c$, a limit cycle will appear and coexist with the SEP.  If the damping is decreased further still, the SEP will cease to exist and only the limit cycle (rotating orbital solution) will remain (corresponding to loss of stability).  Fig. \ref{fig:spdf logS} provides a visual illustration of the different regions of stability as a function of damping parameter.

%The damping $\bar{D}$ also influences the computation domain and computation time significantly. At small $\bar{D}$ the stationary density distribution broadens which demands a larger numerical domain( a larger $x_{2max}$), a longer computation time $t_{max}$. At large $\bar{D}$ the distribution becomes too narrow to numerically resolve it more quickly. We verified that, starting from an arbitrary initial condition, the probability density always converges to the same the stationary density function but requires different computation time.




\section{Correlation to Transient Stability Analysis}



\subsection{Deterministic Stability Analysis}

The stability of the deterministic SMIB model
  \begin{eqnarray}
    \!\!\! \dot{x}_1&=& x_2 \notag \\
     \!\!\! \dot{x}_2 &=&  \bar{P}_m-\sin x_1-\bar{D}x_2   \label{eq:SMIB 3}
  \end{eqnarray}
is briefly reviewed in this section as a precursor to generalization to the stochastic case.
Since (\ref{eq:SMIB 3}) is periodic in $x_1$, the state space $(x_1,x_2)$ can be viewed either in planar space $R^2$ or cylindrical space $S^1 \times R$ where $S^1=[-\pi,\pi]$. The cylindrical space $S^1 \times R$ is a more natural space from the physical perspective. For example, there are only two equilibria in $S^1 \times R$, but there are infinite number of equilibria in $R^2$ (repeated every $2\pi$ radians). Both spaces will be used when convenient. For simplicity, the SMIB parameters are chosen such that there are always two equilibria in $S^1 \times R$: one is the SEP (a sink) with $(\delta_s, 0)$ and the other is a saddle point with $(\delta_{u1},0)$ where $\delta_s=\arcsin(\bar{P}_m)$ and $\delta_{u1}=\pi-\delta_s$.

In this analysis, only the damping parameter $\bar{D}$ will be varied, while $\bar{P}_m$ will remain fixed. Depending on the damping level, deterministic SMIB models can have one of three possible structures for a stability phase plane, as illustrated in Fig. \ref{fig:SMIB_solutions}.
If the damping parameter exceeds the critical damping level (e.g. $\bar{D}>\bar{D}_c$), a typical phase plot is illustrated in Fig. \ref{fig:SMIB_solutions}(a). In this case, the saddle nodes $(\delta_{u1},0)$ and $(\delta_{u1}-2\pi,0)$ define the boundary of the region of attraction of the SEP. The boundaries of the various equilibria are shown as solid lines.  A typical state trajectory is shown with the center dashed-line arrow.  All initial conditions that lie within the region indicated by the middle solid lines will converge to the central SEP.  Initial conditions to the right of this region or above the boundary surface will lead to trajectories that converge to an equilbria that lies on the $\delta$ axis to the right of $2\pi$.  Similarly initial conditions to the left of the central region of attraction will converge to the leftmost  equilibria.  Even though in the planar space $R^2$ there appear to be multiple equilibria, there is only one SEP solution in cylindrical space $S^1 \times R$ due to the periodicity of $\delta$.

The critical case, $\bar{D}=\bar{D}_c$ is illustrated in Fig. \ref{fig:SMIB_solutions}(b), in which a heteroclinic orbit connects the two saddle nodes. As the damping is decreased below the critical damping level, the region of attraction splits and trajectories may be drawn towards either the SEP or the limit cycle based on the initial condition.  Several examples of different trajectories are illustrated by the dashed lines in Figure \ref{fig:SMIB_solutions}(b).  The value of $\bar{D}_{c}$ at which this occurs can be estimated numerically \cite{fp:Alberto00}.


For the situation in which $\bar{D}<\bar{D}_c$, the saddle node equilibria $(\delta_{u1},0)$ sits on the boundary of the region of attraction of the SEP leading to a stable rotating solution (corresponding to a limit cycle on the cylindrical surface $S^1 \times R$). The phase plot shown in Fig. \ref{fig:SMIB_solutions}(c) shows the co-existence of the two attractors: the SEP and the limit cycle.  The dashed lines show several possible trajectories with different initial conditions.




\subsection{Impact of load perturbations on stability}

\begin{figure}[b]
\centering
\includegraphics[height=1.8in]{ode_example.eps}
%\includegraphics[width=1.65in]{sde_example.eps}
\caption{Unstable trajectory in the deterministic system with $\bar{D}=\bar{D}_c=0.409$}
\label{fig:SMIB ODE_Phase}
\end{figure}

\begin{figure}
\centering
\includegraphics[height=1.8in]{sde_example.eps}
\caption{Several individual trajectories in the stochastic system with $\bar{D}=\bar{D}_c=0.409$ }
\label{fig:SMIB SDE_Phase}
\end{figure}

In the deterministic system, systems with damping parameters less the $\bar{D}_c$ have regions of initial conditions for which the SEP cannot be attained.  For example, consider the region $ 0.35<\bar{D}<0.409$ in Fig. \ref{fig:spdf logS}. For a deterministic system, this case corresponds to $ \bar{D}<\bar{D}_{c}$: the SEP and the limit cycle coexist. The initial condition of the system determines to which attractor the trajectory will converge. Fig. \ref{fig:SMIB ODE_Phase} shows that a deterministic trajectory with critical damping $\bar{D}_c$ starting from $x(0)=[0, 2]$ converges into a rotating solution, or limit cycle. This trajectory is unstable in the sense of power system stability.

However, when noise is added to the system, trajectories close to boundaries of the regions can inadvertently cross the boundaries and exhibit behavior not necessarily predicted by their initial conditions in the deterministic system.  This can be traced back to the behavior of the probability density function developed in the last section.  For damping values close to the critical damping value $\bar{D}_c$, the PDF exhibits both ridge and peak characteristics.  As seen from Fig. \ref{fig:spdf logS}, that although both peak and ridge characteristics exist, the peak is prevalent and dominates the probability flow.  Therefore, in a stochastic system, even though the damping may be less than the critical damping, the majority of trajectories will still be attracted by the SEP and not the limit cycle.  In short, the noise in the system actually {\bf improves} the probability of stability.  Fig. \ref{fig:SMIB SDE_Phase} shows a set of the same trajectory as in Fig. \ref{fig:SMIB ODE_Phase} except that a series of Gaussian perturbations have been injected into the system.  Note that while the deterministic system is {\em always} unstable, the stochastic system shows a high probability of attaining the SEP. \hl{It should be mentioned that there is no general conclusion that noise always improves the power system stability. The observation here implies that probability flows near the UEP are subject to crossing the stable or unstable manifold in the presence of system noise. This phenomenon significantly affects the probability attracted by the SEP as the trajectories approach UEP.}



%Improvement of stability is not an infrequent occurrence in stochastic systems.  Consider for example the scalar stochastic process $x_t$ given by the first order Ito stochastic differential equation
%\begin{equation}
%dx_t=rx_tdt+\alpha x_t dW_t
%\label{eq:stoch}
%\end{equation}
%where $r$ is a scalar constant, $dW_t$ is the injected Gaussian noise, and $\alpha$ is the magnitude of the noise.  The explicit solution to this equation is
%\begin{equation}
%x_t=x_0 \exp \left( \left(r-\frac{1}{2}\alpha^2\right)t+\alpha W_t \right)
%\end{equation}
%The qualitative behavior of the process as $t \rightarrow \infty$ is
%\begin{enumerate}
%\item If $r-\alpha^2/2>0$, then $x \rightarrow \infty$ with probability 1.0
%\item If $r-\alpha^2/2<0$, then $x \rightarrow 0$ with probability 1.0
%\item If $r-\alpha^2/2=0$, then $x$ fluctuates between arbitrarily large and arbitrarily small values with probability 1.0
%\end{enumerate}
%
%Note that the stability response is not fully governed by the deterministic boundary $r=0$, but rather that noise may actually improve the stability of the system.  Figure \ref{fig:stoch} shows the solutions to equation (\ref{eq:stoch}) for $r=1$ and values of $\alpha=1, \sqrt{2}, 2$ for the same $dW_t$ in each run.
%
%
%
%\begin{figure}
%	\centering
%		\includegraphics[height=2.0in]{stoch.eps}
%	\caption{Examples of different randomness levels in equation (\ref{eq:stoch})}
%	\label{fig:stoch}
%\end{figure}

%We first summarize the main points of the preceding sections.
%\begin{enumerate}
%  \item Stationary probability density of SMIB model has diverse shapes upon the parameters. It could be peak-shaped, or ridge-shaped, or with more complicated shape.
%  \item Additive noise may destroy the limit cycle attractor and make positive contribution on improvement of the stability under certain conditions.
%\end{enumerate}

%Firstly, the proposed FPE is a two-dimensional parabolic PDE . let us remark that the numerical treatment is not complete.



\section{Application to Multimachine Systems}

\subsection{Principle}
The application of the Fokker-Plank equation to multimachine power systems is extremely challenging due to the high dimensionality of the numerical solution.  However, we propose to extend the Fokker-Plank approach to multimachine systems through the use of the extended transient stability.  The proposed  procedure is shown in Fig. \ref{fig:SME} and outlined below:

\begin{enumerate}
\item {\em Find the multimachine--SMIB equivalent}

Extended transient stability  is an approach by which the generators in the system are separated into two groups which can then be replaced by a two-machine system and solved via a SMIB equivalent \cite{fp:Xue88}-\cite{fp:Pavella94}. This method relies on the observation that the loss of synchronism in a multimachine system is the  result of the machines' irrevocable separation into two groups. Therefore, for a given contingency:
\begin{enumerate}
  \item apply an identification algorithm to determine the machine groupings,
  \item transform the multidimensional multimachine dynamic equations into a two machine equivalent system, and
  %\item assess stability based on equivalent SMIB system: critical clearing time (CCT), or stability margin.
  \item compute the state trajectory and parameters of the SMIB system at any given clearing time (not restricted to the critical clearing time) through time-domain simulation.
\end{enumerate}


\item {\em Calculate the probability of stability of the system}

The dynamic equations of the multimachine system may be transformed into an equivalent SMIB set of equations:
\begin{eqnarray} \label{eq:SMIB equivalent}
M \ddot{\delta}=P_m-P_c-P_{max}\sin(\delta-\nu)-D\dot{\delta}
\end{eqnarray}
where $M,P_m,P_c,P_{max}$ and $\nu$ are parameters which are determined by the contingency and operating conditions \cite{fp:Pavella94}. If any parameters of (\ref{eq:SMIB equivalent}) have a noise component, then the ODE can be reformulated as an SDE. Therefore, the Fokker-Planck equation can be derived to describe the evolution of the probability density function of the critical machine of the multimachine system. The perturbation of $P_m$ and $P_c$ can introduce additive noise term into the system and the perturbation of $P_{max}$ and $\nu$ will introduce multiplicative noise term. \hl{It is noted that}$P_c, P_{max}$ and $\nu$ \hl{ vary with all the bus voltage and angles in the system and hence vary with time. In practice, these time-dependent parameters can be set as "base case" plus "perturbation" value. Perturbation terms can be set as time-dependent functions of other variables exported from deterministic time-domain simulation. Further, it is assumed that these time-dependent parameters are "close" to the base case parameters at the deterministic operating point (i.e., perturbations are relatively small values). In this case, it becomes palatable to set the perturbation as simple constant noise term. Nonetheless, the reduction of large scale systems into SIMB equivalents might cause information loss and error in assessment of stability.}

\vspace*{0.1in}
We propose a \emph{probability of stability} measure $Pr_s$ to represent the probability that the system trajectories remain in a stability region $A$:
\begin{eqnarray} \label{eq:Prs}
    Pr_s(t)&=&{ \iint\limits_{\{ \delta,\omega\} \in A}} p(\delta,\omega,t)d \delta d \omega
\end{eqnarray}

$Pr_s(t)$ is a time-varying variable determined by the transitional PDF $p(\delta,\omega,t)$ and stability region $A$. The stability region is defined as the  region surrounding the SEP (e.g. as approximated by the equipotential energy surface through the controlling UEP). This procedure is summarized:

%Note the definition of stable region will affect significantly the value of $Pr_s$. The stable region can be taken roughly as the rectangle area, or taken as the more elegant stability region surrounding SEP. In this stage,

%$(|\omega|<\omega_{max},|\delta|<\pi)$

%\begin{eqnarray} \label{eq:Pr Stable Def}
%Pr_s(x,t|x_0,0)=Pr\{x: x\in A)\}
%\end{eqnarray}
%If the transition probability density $p(x,t|x_0, 0)$at time $t$ is known, (\ref{eq:Pr Stable Def}) can be written as
%\begin{eqnarray} \label{eq:Pr Stable 1}
%Pr_s(x,t|x_0,0)=\int\limits_{A(x_s)}p(x,t|x_0,0)(dx)
%\end{eqnarray}

%where $\chi$ is the feasibility region of the solution of Fokker-Planck equation. The normalization condition is that $\int\limits_{\mathcal{X}}p(x,t|x_0,0)(dx)=1$.
%Schematically, this stage conforms to the following procedure:
\begin{enumerate}
  \item construct the Fokker-Planck equation using the parameters and initial conditions obtained from the SMIB equivalent system,
  \item solve the corresponding FPE numerically to yield the transitional probability density function, and
  \item assess the evolution of the probability of stability over time $t$ by (\ref{eq:Prs}).
\end{enumerate}

\end{enumerate}


\begin{figure}[t]
	\centering
		\includegraphics[width=3.4in]{SME.eps}
	\caption{Extended transient stability framework based on the Fokker-Planck equation and equivalent  SMIB  }
	\label{fig:SME}
\end{figure}

\subsection{Illustrative Examples}
\subsubsection{Case I (Three-generator WECC System)}
The proposed method is applied to the well-known three-generator WECC test system. The equivalent SMIB is composed of two coherent areas (generators 2 and 3) and (generator 1).  Fig. \ref{fig:CCT} shows the deterministic SMIB equivalent representation for several clearing times. The critical clearing time (CCT) is determined to be $0.17s$ for the given contingency.  Fig. \ref{fig:ProbStab_t}  shows the evolution of the probability of stability index of (\ref{eq:Prs}) for critical clearing time. The trace represents the probability of the system staying in the stability region. It is shown that the probability of stability $Prs$ drops as the probability flows approaching UEP near $1s$ and converge into a stationary value eventually. Probability of stability is a time-varying index and we are more interested in its stationary value. Table \ref{tbl:WECC Results} lists the stationary value of the probability of stability index for clearing times $0.12s, 0.17s$ and $0.2s$. It is noted that the distinction between the stable case and unstable case in the deterministic power system is not clear-cut in the stochastic power system. The proposed stability measure can provide a quantitative measure of stability in the presence of system noise.

 
 %Fig. \ref{fig:ProbStab_t} shows the evolution of the probability of stability index for clearing times  $0.12s, 0.17s$ and $0.2s$.


% The $Pr_s$ traces start at the respective clearing times (not at $t=0$).The trace for a clearing time of $0.12s$ is 1.0 for any time indicating that the system is always stable.  The trace for a clearing time of $0.17s$ is on the boundary of stability.  The probability of stability drops as the state gets captured by the limit cycle shortly after the clearing time.  Note that the probability of stability increases slightly with time because there is the possibility that the noise in the system will push the system state off the limit cycle and towards the SEP since, as noted previously, the presence of noise in the system may improve the stability.  The trace for the clearing time of $0.2s$ shows an initial increase in stability as the generator ``swings in'' immediately after the fault is cleared, but then rapidly decreases again as synchronism is lost.  This example indicates that the proposed stability measure can provide a quantitative measure of stability in the presence of system noise.
%
% The values of probability of stability drop at around $1s$ because the probability flows are closing to UEP. The noise forces some portion of flows out of the stability region of SEP and the make $Pr_s$ less than $100\%$ when time $t$ goes infinity. We only consider the light noise level because the noise level cannot be too much in practical power system. Instead of CCT, the new measure $Pr_s(t)$ provide the quantitative assessment for provability evolution of the states staying in the safe space.

\begin{table}[!h]
\renewcommand{\arraystretch}{1.3}
\caption{Probability of Stability for the WECC system }\label{tbl:WECC Results}
\begin{center}
\begin{tabular}{|l |l | l |}
\hline
Clearing Time& Deterministic     & Probability of  Stability  \\
In Seconds  & Stability & Stationary Value\\
\hline
\hline
\quad \quad 0.12    & \quad Stable  & \quad \quad \quad 0.9999   \\
\quad \quad 0.17 *    & \quad Stable & \quad \quad \quad 0.9628   \\
\quad \quad 0.20    & \quad Unstable & \quad \quad \quad 0.3808  \\
\hline
\end{tabular}
\end{center}
%\newline
%\newline
\begin{tabular}{l}
* deterministic CCT\\
\end{tabular}
\end{table}

%Moveover, the affect of parameters on the probability of stability have been studied. We only consider the light noise level because the noise level cannot be too much in practical power system. Similarly with the deterministic system, the large damping factor can significantly improve the probability of stability. However, different with the deterministic case, the deterministic CCT can not give the precise boundary of the stable and unstable orbits. The new measure $Pr_s$ instead assess probability of the states staying in the stable region.

\begin{figure}
	\centering
		\includegraphics[height=1.8in]{SME_TC.eps}
	\caption{Phase plane: the equivalent SMIB of the three machine system}
	\label{fig:CCT}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[height=1.8in]{ProbStab_t.eps}
	\caption{Probability of stability time evolution}
	\label{fig:ProbStab_t}
\end{figure}


%\begin{figure}
%	\centering
%		\includegraphics[height=1.6in]{ProbStab_D.eps}
%	\caption{Probability of Stability for CCT}
%	\label{fig:ProbStab_D}
%\end{figure}

\subsubsection{Case II (IEEE 50 generator system)}
To further validate the proposed method, the simulations have been conducted on a larger testing system of the IEEE 50-generator, 145-bus system \cite{fp:TSReport92}. The contingencies are specified as three phase fault at Bus 7 cleared by opening line between Bus 7-Bus 6. Fig. \ref{fig:50m} shows the deterministic dynamic trajectories of stable case and unstable case respectively. In the unstable case,there are two generators at Bus 104 and Bus 111 are identified as the critical machines (CMs) which are the group losing the synchronism with other machines. Applying the proposed method, the results of probability of stability can be obtained as shown in Table. \ref{tbl:50m Results}. The critical clearing time can be determined as a time between $0.11s$ and $0.12s$ in the deterministic power system. However, similarly with the previous analysis, the difference between the values of probability of stability for these two clearing time are not rigid and clear. 

\begin{figure}
	\centering
		\includegraphics[width=1.72in]{50m-Stable.eps}
        \includegraphics[width=1.72in]{50m-Unstable.eps}
        \begin{scriptsize}
        \begin{tabular}{cc}
        (a) Stable case, & \hspace{1in} (b) Unstable case
        \end{tabular}
        \end{scriptsize}
	\caption{Phase plane: the equivalent SMIB of the three machine system}
	\label{fig:50m}
\end{figure}


\begin{table}[!h]
\renewcommand{\arraystretch}{1.3}
\caption{Probability of Stability for the IEEE 50-generator system }\label{tbl:50m Results}
\begin{center}
\begin{tabular}{|l |l | l |}
\hline
Clearing Time& Deterministic     & Probability of  Stability  \\
In Seconds  & Stability & Stationary Value\\
\hline
\hline
\quad \quad 0.11     & \quad Stable & \quad \quad \quad 0.9628   \\
\quad \quad 0.12    & \quad Unstable & \quad \quad \quad 0.9510  \\
\hline
\end{tabular}
\end{center}
\end{table}


%\subsubsection{Non-Gaussian Perturbation}

%\hl{Gaussian white noise is natural choice for modeling many random variables which will be normally distributed, such as random load level in power system. However, it should be mentioned that other types of random fluctuations are possible, such as wind energy, short term fluctuations in load/generation. [] proposed to
%
%Typically, SDEs incorporate white noise which can be thought of as the derivative of Brownian motion (or the Wiener process); however, it should be mentioned that other types of random fluctuations are possible, such as jump processes.
%
%
% When the Gaussian white noise is not suitable, the color noise or even more complex random process can be incorporated using the white noise (or Wiener process) as the fundamental process. For example, the Rayleigh process can be modeled as one-dimensional SDE which has the stationary PDF of Rayleigh distribution. And Rayleigh distribution is a particular Weibull distribution that is well-known to model wind energy.
%}

\section{Conclusions and Future Work}
This paper analyzed the stochastic SMIB model under white noise perturbations. The Fokker-Planck equation was used to model the probability density function of the stochastic SMIB model and describe its evolution over time.  A qualitative assessment indicated that the PDF exhibits several dominant characteristics include peaks and ridges that correspond to different stability responses.  These differences were interpreted with respect to a deterministic system and then generalized to the stochastic system with the same defining parameters.  The numerical results verify that the addition of noise may improve  dynamic stability under certain conditions. Although this paper has focused on the relatively straightforward SMIB, it also lays the framework when combining the equivalent SMIB techniques for further research into the stability behavior of large scale stochastic power systems that exhibit nonlinear stability behavior.

%\appendices
\appendices



%\section{Alternating Direction Implicit Finite Difference Methods}\label{sec:FDTD}
%
%Finite difference methods are standard tools in the solution of partial differential equations. In this method, differential operators for space or time are approximated as finite difference with small steps. In this study, we focus on the two-dimensional nonlinear FPE (\ref{eq:SMIB FPK}) for the SMIB model. We use an Alternating-Direction-Implicit (ADI) algorithm to solve the proposed FPE \cite{fp:LeVeque07}. The ADI finite difference scheme has two advantages: 1) the ADI method splits the numerical procedure into two stages which reduces the dimension of the system in each stage. and 2) the implicit scheme allows a larger time step size to be taken.
%
%The numerical domain is $-\pi< x_1< \pi$, and $ x_{min}<x_2<x_{max}$ where $x_{max}$ and $x_{min}$ depend on the SMIB parameters and initial conditions.
%In a general formulation, an ODE system is used to represent the finite difference system of (\ref{eq:SMIB FPK})
%\begin{eqnarray} \label{eq:FDTD ODE}
%  \frac{dP}{dt}&=& (L_1+L_2)P \\
%  P(0) &=& P^0 \notag
%\end{eqnarray}
%where the operator $L_1$ represents the discretization in the $x_1$-direction, and $L_2$ represents the discretization in the $x_2$-direction. The operator splitting (or {\it alternating direction implicit} (ADI)) methods are suggested for use in solving parabolic PDEs in a sequence of stages \cite{fp:LeVeque07}. Each stage involves only one operator. The ADI method can reduce the size of the ODE (\ref{eq:FDTD ODE})  to be solved at each stage.
%
%Stage 1: a half-step is taken implicitly in the direction of $x_2$ and the problem is treated as a 1D subproblems along the line $x_1(i)= const$. The operator $L_2P=\frac{\partial}{\partial x_2}\left[ -f_2 P+\frac{\sigma^2}{2}\frac{\partial}{\partial x_2}P \right]$. We used the implicit scheme of \cite{fp:Zorzano97} which is unconditionally stable for solving parabolic problems with lower order terms :
%
%\begin{flalign} \label{eq:OSM_S1}
%  &\frac{P^{n+\frac{1}{2}}_{i,j}-P^{n}_{i,j}}{\Delta t}= \frac {\sigma^2} {2} \frac{P^{n+\frac{1}{2}}_{i,j+1}+P^{n+\frac{1}{2}}_{i,j-1}-2P^{n+\frac{1}{2}}_{i,j}}{(\Delta x_2)^2} \notag \\
%  &-\frac{f_{2_{i,j+\frac{1}{2}}}(P_{i,j+1}^{n+\frac{1}{2}}+P_{i,j}^{n+\frac{1}{2}})}{2 \Delta x_2} +\frac{f_{2_{i,j-\frac{1}{2}}}(P_{i,j}^{n+\frac{1}{2}}+P_{i,j-1}^{n+\frac{1}{2}})}{2 \Delta x_2}
%    %\\  P^{n+\frac{1}{2}}(t^n) &=& P^n
%\end{flalign}
%
%Stage 2: a half-step is taken implicitly in the $x_1$-direction and the problem is treated as a 1D subproblem along the line $x_2(j)= const$. When the operator  $L_1P=\frac{\partial}{\partial x_1}\left[ -f_1 P\right]$ is applied, the scheme is written:
%\begin{flalign}\label{eq:OSM_S2}
%  &\frac{P^{n+1}_{i,j}-P^{n+\frac{1}{2}}_{i,j}}{\Delta t}= -x_2\frac{P^{n+1}_{i+1,j}-P^{n+1}_{i-1,j}}{2 \Delta x_1}
%\end{flalign}
%
%(\ref{eq:OSM_S1}) and (\ref{eq:OSM_S2}) are equivalent to
%\begin{eqnarray}
%  \left ( I-\frac{\Delta t}{(\Delta x_2)^2}M_1 -\frac{\Delta t}{2\Delta x_2} M_2\right) P^{n+\frac{1}{2}} &=& P^{n} \\
%  \left ( I-\frac{\Delta t}{2 \Delta x_1}M_3 \right) P^{n+1} &=& P^{n+\frac{1}{2}}
%\end{eqnarray}
%where $I$ is identity matrix and $M_1$, $M_2$ and $M_3$ are tridiagonal matrices derived from (\ref{eq:OSM_S1}) and (\ref{eq:OSM_S2}).


% use section* for acknowledgement
\section*{Acknowledgment}

The authors gratefully acknowledge the financial support of the National Science Foundation under projects EFRI-0835995 and ECCS-1068996.  The authors would like to extend their thanks to Alex Stankovic (Tufts University) and Bernie Lesieutre (U. Wisconsin) for their helpful comments and suggestions on this paper.
%
%The authors would like to thank...
% Can use something like this to put references on a page
% by themselves when using endfloat and the captionsoff option.

%\nocite{*}
%\bibliographystyle{IEEEannot}
%\bibliography{FPK_Bib}


\begin{thebibliography}{1}

%\bibitem{fp:Pai89}
%M. A. Pai, {\em Energy Function Analysis for Power System Stability}, Kluwer Academic Publishers, 1989.
%
%\bibitem{fp:Chiang95}
%H.D. Chiang, C.C. Chu, and G. Cauley, ``Direct stability analysis of electric power systems using energy functions: Theory, applications and perspectives,'' \emph{Proc. of the IEEE}, vol. 38, n. 11, pp. 1497-1529, Nov 1995.

\bibitem{fp:Kundur94}
P. Kundur, {\em Power System Stability and Control}, New York:  McGraw Hill, 1994.

\bibitem{fp:Billintion81}
R. Billintion and P. R. S. Kuruganty, ``Probabilistic assessment of transient stability in a practical multimachine system,'' \emph{IEEE Transactions on  Power Apparatus and Systems,} vol. 100, no. 7, pp. 3634-3641, July 1981.

\bibitem{fp:Timko83}
K. J. Timko, A. Bose, and P. M. Anderson, ``Monte Carlo simulation of power system stability,'' \emph{IEEE Transactions on  Power Apparatus and Systems,} vol. 102, no. 10, pp. 3453-3459, Oct. 1983.

\bibitem{fp:Wu83}
F. F. Wu and Y-K Tsai, ``Probabilistic dynamic security assessment of power systems: Part I - basic model,'' \emph{IEEE Transactions on  Circuits and Systems,} vol. 30, no. 3, pp. 148 - 159, March 1983.

\bibitem{fp:Hockenberry04}
J. R. Hockenberr and B. C. Lesieutre, ``Evaluation of uncertainty in dynamic simulations of power system models: The probabilistic collocation method,''  \emph{ IEEE Transactions on Power Systems}, vol. 19, no. 3, pp. 1483 - 1491, Aug. 2004.

\bibitem{fp:Hiskens06}
I. A. Hiskens, and J. Alseddiqui `` Sensitivity, Approximation, and Uncertainty in Power System Dynamic Simulation,''  \emph{ IEEE Transactions on Power Systems}, vol. 21, no. 4, pp.1808 - 1820,Nov. 2006.

%\bibitem{fp:Pai95}
%M.A. Pai, P.W. Sauer, and B.C. Lesieutre, ``Static and dynamic nonlinear loads and structural stability in power system,''  \emph{ Proc. IEEE}, vol. 83, pp. 1562-1572, Nov. 1995

\bibitem{fp:Risken89}
H. Risken, {\em The Fokker-Planck Equation: Methods of Solutions and Applications}. New York: Springer, 1989.


\bibitem{fp:Paganini99}
F. Paganini and B. C. Lesieutre, ``Generic properties, one-parameter deformations, and the BCU method,'' \emph{ IEEE Trans. Circuits Syst. I}, vol. 46, pp. 760?763, June 1999.

\bibitem{fp:Alberto00}
L. F. C. Alberto and N.G. Bretas, ``Application of Melnikov's method for computing heteroclinic orbits in a classical SMIB power system model,''  \emph{ IEEE Trans. Circuits Syst. I}, vol. 47, no. 7, pp. 1085-1089, Jul. 2000.

\bibitem{fp:Chen05}
X. Chen, W. Zhang and W. Zhang. ``Chaotic and subharmonic oscillations of a nonlinear power system,''  \emph{ IEEE Trans. Circuits Syst. II}, vol. 52, no. 12, pp. 811-815, Dec. 2005



\bibitem{fp:Andronov66}
A. A. Andronov, A. A. Vitt and S. E. Khaikin, {\em  Theory of Oscillators.} New York:  Dover Publications, 1966.

\bibitem{fp:Levi78}
M. Levi, C. Hoppensteadt, and W. L. Miranker, ``Dynamics of the Josephson junction,''  \emph{ Quart. Appl. Math.}, pp. 167-198, July 1978.

\bibitem{fp:Oksendal07}
B.Oksendal, ``Stochastic Differential Equations,'' Berlin: Springer, 2007.

\bibitem{fp:Odun-Ayo12}
T. Odun-Ayo, M. L. Crow,``Structure-preserved power system transient stability using stochastic energy functions ,''  \emph{ IEEE Transactions on Power Systems}, vol. 27, no. 3, pp. 1450-1458, Aug. 2012


\bibitem{fp:Dong12}
 Z.Y.Dong, J.H. Zhao, D. J. Hill, ``Numerical Simulation for Stochastic Transient Stability Assessment,''  \emph{ IEEE Transactions on Power Systems}, published online.

 \bibitem{fp:Wang12}
 K. Wang, M. L. Crow, ``Fokker-Planck Equation Application to Analysis of a Simplified Wind Turbine Model,''  \emph{ North American Power Symposium (NAPS) September, 2012}, accepted

%\bibitem{fp:Salam85}
%F. M. Salam and S. S. Sastry, ``Dynamics of the forced Josephson junction circuit: The region of chaos,'' \emph{ IEEE Trans. Cirrcuits Syst.}, vol. 32, no. 8, pp. 784?796, Aug. 1985.

%\bibitem{fp:Molina-Garcia11}
%A. Molina-Garcia, M. Kessler, J. A. Fuentes, and E. Gomez-Lazaro ``Probabilistic characterization of thermostatically controlled loads to model the impact of demand response programs,'' \emph{ IEEE Trans. Power Systems}, vol. 26, no. 1, pp. 241-252, Feb. 2011.

%\bibitem{fp:Arnold77}
%L. Arnold, {\em Stochastic Differential Equations: Theory and Applications}, John Wiley \& Sons, 1977.



%\bibitem{fp:Feller54}
%W. Feller, ``Diffusion process in one dimension,'' {\em Transaction of American Mathematical Society}, vol. 77, pp 1-31, 1954.

%\bibitem{fp:Fuller69}
%A. T. Fuller, ``Analysis of nonlinear stochastic systems by means of the Fokker-Planck equation,'' {\em International Journal of Control}, vol. 9, pp 603-655, 1969.
%\bibitem{fp:Zhu96}
%W. Z. Zhu and Y. Q. Yang, ``Exact stationary solutions of stochastically excited and dissipated integrable Hamiltonian systems,'' {\em ASME J. Appl. Mech.}, vol. 63, pp. 493?500, 1996.

%\bibitem{fp:Zhu90}
%W.Q. Zhu, G. Q. Cai, and Y. K. Lin, ``On Exact Stationary Solutions of Stochastically Perturbed Hamiltoian Systemss'' {\em Probabilistic Engineering Mechanics,} vol. 5, no. 2, pp 84-87, 1990.
%
%\bibitem{fp:Zhu91}
%W.Q. Zhu, G. Q. Cai, and Y. K. Lin, ``Stochastically Excited Hamiltoian Systemss'' {\em Proceedings of IUTAM Symposium on Nonlinear Stochastic Mechanics,} eds. N Bellomo and F. Casciati, pp 543-552, 1991.

%\bibitem{fp:Arnold98}
%L. Arnold, {\em Random Dynamical Systems}, Springer-Verlag, 1998.
%
\bibitem{fp:Kloeden99}
 P. E. Kloeden, E. Platen, {\em Numerical Solution of Stochastic Differential Equations, } 2nd ed. Berlin:  Springer-Verlag,  1999.

 \bibitem{fp:Higham01}
D. L. Higham, ``An Algorithmic introduction to numerical simulation of stochastic differential equations,'' \emph{SIAM Review}, vol.43, no.3, pp.525-546, 2001.

 \bibitem{fp:Strogatz01}
Steven H. Strogatz, {\em Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering (Studies in Nonlinearity)}, Westview, 2001

\bibitem{fp:Lin04}
Y. K. Lin, G. Q. Cai, {\em Probabilistic Structural Dynamics. } New York:  McGraw-Hill, 2004.


% \bibitem{fp:Zorzano99}
%M.P. Zorzano, H. Mais, L. Vazquez. "Numerical solution of two dimensional Fokker-Planck equations" , {\em Applied Mathematics and Computation }, vol. 98, pp. 109-117 , Feb. 1999.

%\bibitem{fp:Zorzano97}
%M. P. Zorzano, H. Mais, L. Vazquez. ``Numerical solution for Fokker-Planck equations,'' {\em Proceedings of the Particle Accelerator Conference}, May 1997.

%\bibitem{fp:Richtmyer67}
%R. D. Richtmyer, K. W. Morton. {\em Difference Methods for Initial-Value Problems}, Interscience Publishers, 2nd ed. 1967.

\bibitem{fp:LeVeque07}
R. J. LeVeque, {\em Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems}, Philadelphia, PA:  SIAM, 2007.

\bibitem{fp:Xue88}
Y. Xue, T. Van Cutsem and M. Pavella, ``A simple direct method for fast transient stability assessment of large power systems,'' \emph{ IEEE Trans. Power Systems}, vol. 3, no. 2, pp. 400-412, May. 1988.

\bibitem{fp:Pavella94}
M. Pavella, P.G. Murthy {\em Transient Stability of Power Systems: Theory and Practice}, UK: Wiley, 1994.

\bibitem{fp:TSReport92}
IEEE Stability Test Systems Task Force of the Dynamic System Performance Subcommittee, `` Transient stability test systems for direct stability methods,'' \emph{ IEEE Trans. Power Systems}, vol. 7, no. 1, pp. 37-43, Feb. 1992





\end{thebibliography}


%\begin{IEEEbiography}
%{\bf Keyou Wang} (M'09) received his B.S. and M.S. degrees in Electrical Engineering from Shanghai Jiaotong University, China in 2001 and 2004 respectively, and his Ph.D. degree from the Missouri University of Science \& Technology (formerly University of Missouri-Rolla) in 2008.
%
%He is currently a Research Fellow of Electrical Engineering at the Missouri University of Science \& Technology. His research interests include power system dynamic and stability, flexible ac transmission system (FACTS) and renewable energy integration.
%\end{IEEEbiography}
%
%\begin{IEEEbiography}
%{\bf M. L. Crow} (S'83--M'90--SM'94--F'10)
%received the B.S.E. degree from the University of Michigan, Ann Arbor, and the
%Ph.D. degree from the University of Illinois, Urbana/Champaign.
%
%She is currently the
%Director of the Energy Research and Development Center and the F. Finley
%Distinguished Professor of Electrical Engineering at the Missouri University of Science \& Technology.  Her research interests include  computational methods for dynamic
%security assessment and the application of power electronics in
%bulk power systems.
%\end{IEEEbiography}

\begin{IEEEbiographynophoto}
{\bf Keyou Wang} (M'09) received his B.S. and M.S. degrees in Electrical Engineering from Shanghai Jiaotong University, China in 2001 and 2004 respectively, and his Ph.D. degree from the Missouri University of Science \& Technology (formerly University of Missouri-Rolla) in 2008.

He is currently a Research Fellow of Electrical Engineering at the Missouri University of Science \& Technology. His research interests include power system dynamics and stability, flexible ac transmission system (FACTS).
\end{IEEEbiographynophoto}

\begin{IEEEbiographynophoto}
{\bf M. L. Crow} (S'83--M'90--SM'94--F'10)
received the B.S.E. degree from the University of Michigan, Ann Arbor, and the
Ph.D. degree from the University of Illinois, Urbana/Champaign.

She is currently the
Director of the Energy Research and Development Center and the F. Finley
Distinguished Professor of Electrical Engineering at the Missouri University of Science \& Technology.  Her research interests include  computational methods for dynamic
security assessment and the application of power electronics in
bulk power systems.
\end{IEEEbiographynophoto}

\end{document}


